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ABSTRACT 





A triquadratic isoparametric solid element is developed to study the behavior 
of thick isotropic and laminated composite plates. The element is a 27 noded La- 
grangian element based on three dimensional elasticity. Material characteristics are 
accounted by either using laminate plate theory or three-dimensional anisotropic 
theory. Element matrices for nonlinear stability analyses are derived based on total 


Lagrangian formulation. 


Results are presented to compare with analytical solutions to validate the 
elements behavior. The effects of various integration schemes on the element per- 
formance are presented. Convergence studies for laminated composites for different 
fiber orientations are provided to illustrate applications. An analysis of thin plates 
is carried out and results for thick plates are compared with available higher order 


plate theories. One row of elements in the thickness directions gives satisfactory 


results for thick laminates. 
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I. INTRODUCTION 


A. OVERVIEW 

The finite element method provides a general tool to solve problems of contin 1a 
such as heat conduction and fluid flow, but it is most widely used in structural 
mechanics. In structural mechanics, the methodology is applicable for static and 
dynamic response of structures and in predicting the elastic stability limits. 

The focus of the present study is to develop tools to analyze thick laminateu 
composite plates and validate the model >y comparing with known solutions. More 
specifically, the objective of the present study is to develop a finite element for both 
linear and nonlinear analysis using three dimensional elasticity relations. . 

By adopting such theory for thick plates, both isotropic and composite, the 
solutions account for transverse shear stresses. This approach eliminates the limi- 
tations imposed by classical plate theory based on Kirchoff-Love hypothesis (Batoz, 
1950] or higher order shear deformation theories (Reddy, 1984, Lo et al., 1977]. 


B. LITERATURE REVIEW 


In this section, some literature pertaining to the analysis of thick composite 
plates is reviewed. The finite element method has been increasingly used as a 
research tool, as well as a design analysis tool, and the methodology is rapidly 
evolving along with the development of faster and more efficient computers. Basic 
concepts of the theory of finite element analysis are well documented (Cook, et iu., 
1989]. Yang (1987) describes various two dimensional higher order elements as well 
as three dimensional solid elements. Bathe (1982) discusses the general formulation 


of finite elements in nonlinear anaiysis for one, two and three dimensional elements. 


based on the total Lagrangian formulation and the principle of virtual displacements. 
A good source for continuum formulation may be found in Malvern (1969). 

Tsai and Pagano (1968) establish a notation in which composite lamina prop- 
erties are invariant with respect to lamina direction. The laminate theory is well 
documented by Vinson (1987), where the elasticity solution for “structures com- 
posed of composite materials” is given for various cases, such as bending of thin 
plates. Based on laminate theory, Hoskin, et al. (1986) outline procedures involved 
in manufacturing composite components and presents some of its applications. 

A higher order shear deformation theory of laminated composite plates was 
developed by Lo, et al. (1977). A higher order nonlinear theory of thick plates was 
suggested by Reddy (1984a, 1984b, and 1985) and presented solutions (Reddy, 1987) 
and compared numerical results to Pagano’s (1969) elasticity solution for the case of 
cylindrical bending. Other elasticity solutions are given by Timoshenko (1951 and 
1959) and Eisley (1989) who discusses the elasticity solutions. The Heterosis finite 
element was suggested by Hughes, et al. (1978) for thick and thin plate bending 
problems. The effect of reduced integration in isoparametric finite elements was 
presented by Zienkiewicz, et al. (1971). 

In recent years, much work is concentrated on the analysis of buckling and post- 
buckling response of laminated plates and shells using nonlinear analysis. Ramm 
(1982) applies degenerate finite elements to solve buckling of thin shells. Arnold, 
et al. (1983) presents a theoretical analysis procedure for prediction of buckling 
and post-buckling in laminated composite plates and compares the results to exper- 
imental results. A combined numerical and experimental study of the post-buckling 
behavior of composite panel is performed by Natsiavas, et al. (1987). Gujbir et 
al. (1989) use an eight noded biquadratic element to study the effects of transverse 


shear on the stability of laminated plates. Some solution algorithms for nonlinear 











analysis of structures by adapting modified Newton-Raphson and arc-length meth- 
ods are given by Kolar et al (1985) and Ford et al (1987). In the literature reviewed, 
there appears to be no discussion on the higher-order solid element for the analysis 
of thick laminated plates. 

This research addresses the problem of using a tri-quadratic displacement field 
based finite element based on three-dimensional elasticity equations. A total La- 
grangian formulation is used to derive relevant element nonlinear matrices, and 
numerical examples are included to validate the linear portion of the development. 

Analysis of typical examples include slender bars under traction and bending 
loads, thin and thick plates under bending loads and effects of various integration 


schemes. 


C. THESIS OUTLINE 

This section provides an overview of various chapters of the thesis. The total 
Lagrangian formulation for analyzing structures composed of three-dimensional el- 
ements is presented in Chapter II. Element matrices are derived for both linear and 
nonlinear static analysis using the incremental load method. The material charac- 
teristics account for both linear isotropic and anisotropic behavior. Formulas are 
provided to obtain work-equivalent loads for distributed body and surface forces. 

Chapter III addresses aspects of computational implementation of the problem 
formulated in Chapter II. Test cases, example calculations and comparison with 
classical solutions and other high order theories are given in Chapter IV. Finally, 


Chapter V summarizes the results and reflects some suggestions for future work. 

















II. THEORETICAL FORMULATION 


A. INTRODUCTION 

In this chapter, using the principle of virtual displacements, the stiffness matrix 
will be developed for static equilibrium of triquadratic isoparametric solid elements. 
In the total formulation presented, both small and large displacements are permis- 
sible for linear and nonlinear structural analysis. For both cases, small strains and 
linearly elastic material will be assumed. 

The element is developed for analysis of both isotropic and composite struc- 


tures. 


B. GENERAL DERIVATION OF FINITE ELEMENT EQUILIBRIUM 
EQUATIONS 
The principle of virtual work is invoked for the general formulation of equilib- 
rium [Bathe, 1982; Cook, 1989]. The principle of virtual work states that a body is 
in equilibrium, if and only if, the total virtual work done by the internal forces is 


equal to the total virtual work done by the external forces. That is, 
bWint = OWest (2.1) 


This principle is equivalent to the minimum total potential energy principle 
(611, = 0}, and holds at any given time. 

Consider a three-dimensional body under arbitrary loads as shown in Figure 
2.1. Using a Cartesian system, let the loads be given by 


Taal (2.2) 














Vid ane ae (2.3) 


(F}} = [FF Fl’ (2.4) 
where {f*}, {f°} and {F"} are surface tractions, body forces, and concentrated 
applied forces respectively. 

The displacements of a finite element in the body due to external load is 


denoted by {d}, where 
{d} =[uv wv)? (2.5) 


and the corresponding strains are given by, 

{€} = [ecs ty Ess Cys Exe Exy |” (2.6) 
for which the corresponding stresses are, 

{0} = [ore Oyy Ose Sys Ses Fey]? (2.7) 


The total internal virtual work for a finite element in the body is {5¢}7{o}dv and 
for the whole body, 
5Wine = | {6e}" {o}dv (2.8) 


where the virtual strains, {5¢}, are 
{5} = [S€ce Sey S€s2 5€ye Stes SExy)” (2.9) 
The total external virtual work is given by: 


&Weee = | {6d}7 {f3}dv + i {5d°}" {f*}ds + {60°} {F} (2.10) 
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Figure 2.1: General 3-D Body 








where {5d’} denotes a surface displacements and {éd"} represents point (nodal) 
displacements corresponding to the applied loads, and the virtual displacements 
{6d} are 

{6d}? = [6u bv dw] (2.11) 


On substituting equations 2.8 and 2.10 into 2.1, we get, 
[{5e}? {o}av = [i5a}? (f8}av + [ {5a"}7 (phd + [6d }T{F} (2.12) 


It may be noted that the principle of complementary virtual work may have been 
undertaken, assuming small virtual stresses with true displacements, yielding in an 
analogous expression for equation 2.12. 


Introducing the generalized Hooke’s law for material constitutive relations, 
{o} = [E] {e} + {20} (2.13) 


where {o,} denotes the element initial stresses and [£] denotes the elasticity matrix 
of the element material. 


In general, the strain-displacement relations are given by 


{e} = [B] {d} (2.14) 
while the virtual strains are given by 


{5e} = (B] {6d} (2.15) 


Substituting equations 2.13 and 2.15 into equation 2.12 and simplifying, we have, 


[[{5ay? ((BI" ETB) {a}dv =f {5a)7(fP}dv 4 f {5ar}* (Fhe 


| {Sd}? [BI] {o,}du + {6d°}7{F"} (2.16) 


The integrations in equation 2.16 are performed over the element volume and 


surface, i.e., we can evaluate every integral using the element local coordinates and 


7 








assemble tor the global system coordinates. Thus, we define the global displacement 


vector and the global virtual displacement vector as follows: 
{D} = [uw uv w, 1s UnUnWn] (2.17) 


and 
{5D} = [6u,5v,6w, ...5u,dv,dwa)” (2.18) 


where n is the total number of nodal points in the body. Now we define, for m 


elements, 
(K] = ¥ [te (eleld» (2.19) 
[e] = [BI (e)[B}av (2.20) 


where [K] and [k;] are the global and local stiffness matrices respectively. In addi- 


tion, we define, 


(Roy = ¥ [IM (s?}dv (2.21) 
{ra}; = [mt (1?) (2.22) 
{R.} = 3 fie? (rhs (2.23) 
{ra}5 = [INT (fds (2.24) 
(Ri) = ¥ [BF {o}do (2.25) 
: {rj = [BF toate (2.26) 


where {R} and {r} denote the global and local load vectors and [N] and (N’] are 
the displacement interpolation (shape functions) matrices for the volume and surface 


where traction is prescribed. Using these definitions, we obtain 


{6D}" [K]} [D] 
{R} 


{6D} ({Ra} + {Ra} - {Rr} + {F'}) (2.27) 
{Rs} + {R,} -— {Rr} + {F"} (2.28) 


8 








By invoking the principle of virtual displacements and noting that {6D} is 


arbitrary, we get the equilibrium equations in the following form: 
[K] {D} = {R} (2.29) 


Equation 2.29 is the basic equation for static equilibrium, which also gives the 


general form for nonlinear analysis with large displacements and strains. 


C. INTERPOLATION SCHEME 
1. Shape Functions (Displacement Interpolation Functions) 
In this section, the interpolation scheme for a triquadratic isoparametric 
solid element will be developed. 
The one-dimensional Lagrange interpolation function based on parame- 
ters is given by 


qg 
P= > MP, = NP, + N2P2+...+ NP, (2.30) 


izl 


where N,, also called the shape functions, are given by 


(2.31) 


A triquadratic solid element is a three dimensional element in which the 
displacements u, v, and w are interpolated by quadratic langrangian interpolation 
functions with 27 nodes. Figure 2.2 depicts an element in the local non-dimensional 
coordinates (r,s, t). 


For an isoparametric element, the geometry may be interpolated as, 


27 
r= >: N;z; 
eS 
y=). Ny 
iz 


9 
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Figure 2.2: Lagrangian Solid Element - 27 nodes 
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27 
z=) Nz; (2.32) 


1=1 


or, in a matrix form, 


Zz 
| y - [N] [ziyizy -.-zeryarzar)” (2.33) 


Zz 


where the shape function matrix is given by 


Nn 0 0 Nz 0 0 ... Nx 0 0 
[N] = 0 N, 0 0 N, 0 eee 0 No 0 (2.34) 
0 0M 0 0 NM... 0 DO Ny 





The shape functions and their derivatives in local coordinates are pre- 
sented in Appendix A. 
2. Jacobian Transformation Matrix . 
To obtain equilibrium equations in the global coordinate system and con- 
struct stiffness matrices, we need the derivatives of the shape functions in cartesian 


system. Using the chain rule for differentiation, we obtain, 





Ni» | Ze Yr 2 Niz 
Nin }=| Ze Yo 2 Niy (2.35) 
Nit Te Ye 2 Ni 
where [J], the Jacobian matrix is given by, 
Ze Vr Zr 
[J] = | Ze Ve. Be | (2.36) 
° Ze Ye 2 


A comma denotes differentiation, where for example, Ni, = oN etc. Using the 
shape function derivatives, the elements of the Jacobian matrix may be calculated 


and is given in Appendix B. The global cartesian derivatives may now be obtained 


Nis Ni 
Niy = | Nino 
Niz Nis 
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as, 


(2.37) 














where the inverse of the Jacobian, 


(T] = (yy (2.38) 
is given explicitly in Appendix B. 


D. STRAIN DISPLACEMENT RELATIONS - [8] 
1. Basic Formulation 
In this section, the basic formulation for nonlinear analysis of a general 
solid body is presented (Refs. Bathe (1982), and Malvern (1969)]. First, some defi- 
nitions and notations will be introduced concerning the coordinate system, displace- 


ment, stress and strain measures and later on the linearized equilibrium equation 


will be developed based on section II/B. 


Consider the motion of a body, or an element within, in a fixed cartesian 
coordinate system as shown in Figure 2.3. We have the body at time 0,t and t + At 
for which the upper left superscript corresponds. The displacements at time t and 
t + At are given as 


Sus =' zy; -° zy (2.39) 


tHAty t+dt 7 0 2. (2.40) 
so that the incremental displacements are 
Uy =tt+4t U; -' U; (2.41) 


where, 


U=uU U=v UuU=Ww 


Ty=r %4=38 t=t (2.42) 


12 








we use the following notation for derivatives at time, say t + At, with respect to 


coordinate at time 0 as, 
gttat,. 


Mig = Or; 
j 


In the present approach, we use the total Lagrangian formulation, refer- 


t+ At 
0 





(2.43) 


encing all variables to the undeformed configuration at time 0, [Ref. Bathe, 1982; 
Malvern, 1969]. It is assumed that at time 0 and ¢ the equilibrium configuration is 
known. Basically, equation 2.12 needs to be solved corresponding to time t + At. 
Since we assume large displacements, and nonlinear constitutive relations [equa- 
tion 2.14], equation 2.12 may be solved by incremental load methods (Ref. Ford & 
Stiemer, 1987]. 

On introducing the 2nd Piola-Kirchhoff stress, it may be shown that the 
2nd Piola-Kirchhoff stress tensor is energentically conjugate to the Green-Lagrange 


strain tensor (Ref. Bathe, 1982; Malvern, 1969]. 
i {et Ate}T {OS }dv = fia.,£8e¥" {tata }ttatdy (2.44) 
where the 2nd Piola-Kirchhoff stress at time t is defined as, 
to} = [tx] {5} [fx] det bx) (2.45) 


such that (‘x], the deformation-gradient tensor is a tranformation operator from the 
coordinates at time 0 to time t. Note that in the equation 2.44, the right hand side 
represents internal virtual work at time t + At over the volume at that time while 
the left hand side has the virtual work integrated over known configuration at the 
reference volume. 

Assuming linear material behaviour, we may use the linear stress-strain 


relations (Generalized Hooke’s Law) for the 2nd Piola-Kirchhoff stress tensor. 
(‘s] = (21 [‘e] (2.46) 
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Figure 2.3: Motion of a body in a fixed Cartesian coordinate system 
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where the Green-Lagrange strain tensor at time t is defined as, 


‘ej = ; (tu; +! uy +! ue, tuk.s) (2.47) 
with i, j, k = 1, 2, and 3. On subs ..uting ‘+4¢u =' u; + uj, we obtain at time 
t+ At, 

Ar ACS. ; (‘u;,. t+ uj; +! us, ‘uns) 

+ ; (ui, + Uj; + ux Ung + Uk, ‘4es) 

+ tei Uk,j (2.48) 
which may be written as, 

Hates a eis + Ei; (2.49) 


where, ‘e,; is defined earlier and the incremental strain ¢;; is given by 
G5 = ei + 5 (2.50) 
In matrix notation, 
{e} = {e} + {n} (2.51) 


The linear incremental strain is identified as 
1 
ej = 2 (ui + ujy +! Ue Uey + UR, ‘uks) (2.52) 


in which ‘u,, and ‘u,; are the known displacement gradients at time t. The non- 
linear incremental strains, then, are givea by 


1 
ny = ati Uk (2.53) 


Rewriting the equilibrium equation as stated in equation 2.12, using the total La- 


grangian approach, we have, 
Joa {Se}T {(+4to}ttOtdy attOt R (2.54) 
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where ‘+°'R, the external virtual work, is assumed to be deformation independent. 
Using the identity form equation 2.44, we may write the equilibrium equation in the 


undeformed configuration as 
[ {erate} T{ats}dv =" R (2.55) 
Noting that {‘e} is displacement invariant, 
{5'+4e} = {5e} (2.56) 
The 2nd Piola-Kirchhoff stress at time t + At may be expressed as 
{tts} = {'s} + {5} (2.57) 


where {s} is the incremental 2nd Piola-Kirchhoff stress. On substituting Equations 
2.56, 2.51, and 2.57 into 2.55 yields, 


i {Se} {5}°do + [ _{6e}7{'s}°dv + [ {én} ?{'s}°dv (2.58) 
The incremental stresses are expressed using equations 2.47 and 2.57 as 
{s} = [E] {te} — [E] {‘e} (2.59) 
which in view of Equation 2.50 yields 
{s} = [E] {e} (2.60) 
Referring to Equation 2.51 and neglecting the nonlinear strain contribution, we get 
the linearized approximation as 
{s} = [E] {e} (2.61) 


and 


{de}? = {5e}7 (2.62) 
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On substituting these into Equation 2.58 and rearranging yields the linearized in- 


cremental equilibrium equation, 


I _{6e}? [El {e}°du + [ {én}? {'s}°du =" B— [ {6e}7{'s}°dv (2.63) 


2. General Nonlinear Discretization 
The general nonlinear finite element discretization for the 27-noded ele- 
ment is presented based on the total Lagrangian formulation discussed in the pre- 
vious section. Equation 2.52 for the linear incremental strain in cartesian form 


yields 
t t 
Cer = Uztiuz,upt v,u,t' w, we 
e, = uytiu,u,t+'v,v,+'w,w 
vv a VW vw yy wu VY WV WY 
e,, = u,t'u,u,tiv,v,tiw,w 
2 2 z Us 2 Uz 2 Ws 
t t t 
Qe, = vstuytiu,u, tu, ‘ust vyu,tuy ‘vt iw, wstwy ‘ws 
t t 
ere = Us t Wet Us te tte Us + Ve Vs te Vet West We 'W, 
t t t t t 
Zesy = Uy tet us Uy tus Uy + Vs Vy ts Vy + We Wy + wes Wy 


(2.64) 
which in matrix form is given by 


{e} = {exo} + {ex} (2.65) 


. 


The first term on the right hand side is displacement independent while the second 


term is displacement dependent with the engineering strains {¢} represented by 
{e} = [ess Cyy Crs 2eys 2€z2 ery]? (2.66) 
We define the incremental displacement gradient in the global coordinates by 
{Ug} = [us Uy Us Ve Vy Vis We Wy ws]? (2.67) 
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Equations 2.64, 2.65, 2.66, and 2.67 result in 


{exo} = [Azo] {Uc} 


{ex} = [Axi] {Uc} 


so that [Azo] and [Az] are given by, 


1 0 
0 0 
0 0 
[Axo] = 
0 0 
0 0 
0 1 
‘tu. O 0 
0 ‘uy, 0 
0 O ‘u, 


‘us, O ‘us 
‘uy ‘us 0 


oo 


oo ° 
bt 
oo 


t t 

Us Uy 
t 

0 ‘uv, 


tv. 0 


oo © 
Oo 


Wy 


0 
0 
1 
0 
0 
0 
0 0 
‘w, 0 
0 ‘w, 
‘we. ‘wy 
0 ‘tw, 
‘we. 0 


(2.68) 


(2.69) 


(2.70) 


(2.71) 


It may be noted that the values of ‘u; ; are known at the new configuration at time 


t+ At. With the displacements interpolated by 


uU 


w 


27 

» N RUE 
kml 

27 

> Nevp 
k=l 


27 
Si: Nywe 
keel 


the local displacement gradients are obtained from 


Ur 


Us 


Us 


a7 
> New Ue 
kw 


27 
> Nie UR 


kul 


27 
> New UR 


k=l 
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(2.72) 


(2.73) 











and similar expressions may be attributed to y and w. These gradients may be 


represented as 


1 


Nie (0) 0 Na+ 0 0 No7° (0) it) wy 
Ur Ni, 0 () Nas 0 0 Na7,0 0 0 
ue Nie (t) 0 No.s o 0 No7,1 0 i) uz 
ue “2 
Ur {0} Nive 0 0 Nar (0) 0 N37," ct) un 
Vs = C9] Ny.4 0 (9) N2,. 0 0 Na7,. 0 
Vie 0 Nye 0 0 Nae 0 0) Nar. 0 
we 
Wie 0 oO Nir ct] 0 Nar 0 tH) No7,° 
we 0 C0] Nis 0 c?] N2,, it) tv) N27, 

0 oO Nis 0 OQ Nas i) 0 N7,6 war 
v27 
war 
(2.74) 

or alternatively, 
{Ux} = [DH] {a} (2.75) 
where the nodal displacement vector is given by, 
{d} = . 2.76 
} = (ur v, wy ug v2 +++ w7] (2.76) 
and the local incremental displacements gradient are given by 
{U,} = z 2.77 
L} = [Up Us Ue Vp Vz Ve Wy Wi, Wy] (2.77) 


As previously mentioned, in isoparametric finite elements, the same inter- 
polation functions are used for approximating the geometry and the displacements. 
Using these definitions, we may transform the displacement in global and local co- 
ordinates by similar transformations used for the geometry. Furthermore, we can 
define arbitrarily the global and local coordinates to coincide at time 0 configuration, 
thereby, the transformation from local coordinates at time 0 to local coordinates at 
any other time, say t or t + At is identical to the transformation that relates the 


global coordinates to the local coordinates at any configuration. In other words, we 
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can use the same jacobian matrix defined previously for all configurations. Writing 


the relation in accordance to equations 2.37 and 2.38 we have, 


{Uc} = [T's] {uz} (2.78) 
such that, 
(T] {0} (0 
[fs] =} (0) [F} (0) (2.79) 
(0) fo) {r) 


and substituting equation 2.75 into 2.78 yields, 


{Uc} = [I's] [DH] {4} (2.80) 


The incremental strains, then, may be expressed in terms of nodal dis- 


placements, and substituting equation 2.77 in 2.68 and 2.69 
{ezo} = [Aco] [T's] [DH] {4} (2.81) 


{ez} = [Ars] U's] [DA] {4} (2.82) 


The strain displacement operator may be identified as 


[Bro] = [Azo] [T's] [DH] (2.83) 
[Bra] = [Axi] [T's] [DH] (2.84) 
such that, 
[Bx] = [Bro + (Bri) (2.85) 
and, 
{e} = [Bz] {4} (2.86) 


It should be underscored that using only the displacement independent strains in 


equations 2.63 and 2.64 results in the linearized problem (same as linear small 
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displacement - small strain formulation) with 


Mz 9 0 Nz, 0 0 Nars 
0 My 0 0 Nry 0 0 
0 0 M. 0 0 Naz, sone 0 
[Bro] = 
0 M,: My 0 Naz Noy 0 
NM, O Me Nz, 0 No N72 
My Mz 0 Nyy Naz 0 Nory 


The displacement dependent contribution to the strain-displacement operator 


is given by 
‘ws Mas "ve Ma ws Me ‘ws Ms 
‘uy M, ‘uy My Wy, Ma ‘uy Ny, 
_ ‘wa Me ‘os Mas wi Ms "ws Nan 
[Bui] = 


wy Matus Me 
"we Ma tus Mo 
‘w, M, +' By Ma 


'y Mate, Mi, 
‘og Mat, Mao 
'v, Ma +' Vy Ma 


‘we Matte. M, 
"ws Mia +' we, Nie 
‘wn Ma + wy Me 


‘wy Ny, +'u, No, 
's, Ny, +‘, Ny 
‘ws Ma, +' By Mas 


To obtain the incremental nonlinear strain contribution, 


tion 2.53, which has the cartesian components, 
1242 2 
nee = 5(u +02 + u4) 
tw = 5(u tut +h) 
‘ss = s(t +v2+w>) 
lys = 5 (tats +u,u, + wyw:) 
Nee = (Uets +20, + we.) 


2 


Ny = g(Uaty + usvy + W2Wy) 
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0 0 
Nary 0 
0 Nary 
No.2 Novy 
0 No72 
N27,2 0 
(2.87) 





‘wa Nas 
‘ws Nig 
‘eo, Nas 

° 
'w, Nyx, +'w, Navy 
‘ws Naw t' 0, Neve 
‘we, Nag t' 0, Nave 


(2.88) 


consider equa- 


(2.89) 


With the variations at time t given by, 


ones = buy: ‘ul + bv,. ‘ve + bw: ‘we 


iny = Sury ‘uy t+ dv, ‘v, + dw, ‘wy 

dn, = buy, ‘u,, + 6v,, ‘v,, + bw, ‘w, 
Qénys = Suy‘u, + ‘uydu, + dv, ‘v, + ‘v dv, + dw, ‘wrt ‘wow, 
26n,, = bu,‘'u, + ‘updu, + dv, ‘v, + ‘v6u, + dw, ‘w+ ‘w rdw, 


Qiney = duz'uy t+ ‘uzdu, + dv, ‘'v, + ‘vsdv,y + dw. ‘wy + ‘w.6dw(2.90) 


(2.91) 


It is worth noting that equations 2.90 are in exactly the same form as the displace- 
ment dependent strains {€z;} given in equations 2.64 and 2.65 with the incremental 


displacements derivatives replaced by their variations. Thus, we may write, 
{6n} = [Ax] {65Ue} (2.92) 
{n} = [Aus] {Uc} (2.93) 
such that the nonlinear incremental strain variation vector is defined as, 
{5n}7 = [5nes Snyy Snes 25nys 26nes 25Nzy] (2.94) 
and [Az)] is as given in 2.71. Observing relation 2.80 for {Uc}, we have, 
{5Uc} = [P's] [DH] {6d} (2.95) 


and substituting for the global variations into the nonlinear strains in equation 2.90, 
we get 


{én}? = {6d}? [DH] [Ts] [Aus]” (2.96) 
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Using the strain-displacement relations defined thus far, we formulate the 
linearized incremental equilibrium equation, given by 2.63 to take the general form 
as stated in 2.29, to give, 

( {tec} ra [*kav|) {d}@+) gs {*4'R} ws! {it py) (2.97) 
where [*kz], in view of equation 2.85 is seen to be 
[tes] = i) [Bx]* [E] [Bz]? dv (2.98) 
which is the linear stiffness matrix, and 7 is the iteration number. This includes the 
displacement dependent and independent contributions. 

When small displacements are assumed, the ['k,] reduces to standard lin- 
ear stiffness matrix, given by fo,[Bzo]” [E]{Bzo]’ dv. In what follows, the derivation 
of [kw z] is described. 

The 2nd-Piola stress vector at time t is accumulated such that, 

{'s} = {'*9*s} + [E] {e} (2.99) 
Using equations 2.51, 2.85, 2.86, and 2.91, the incremental strains take the form 
{e} = (2[Azi] + [Azo]) [T's] [DA] {4} (2.100) 
or alternatively, 
{e} = (2[Bri] + [Bro]) {4} (2.101) 
Noting that equations 2.99 and 2.100 are valid at any time t, and by using equation 
2.46, the 2nd Piola-Kirchhoff stress may be written as 
{'s} = [E] (2[Axi] + [Azo}) [T's] [D4] {4} (2.102) 
and the nonlinear part of equation 2.63 becomes, using the relations 2.95 and 2.101, 
[ony t's}eav =f ({8d}7 [DAY (P]” [Aual”) [2] 
(2[Ari] + [Azol) Ta] [DH] {d}°dv (2.103) 
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we define the nonlinear strain contribution, [kwz], as 
fewe] =f (DAY Ws}? (Anal” (£)(2(4na) + [4rol) [Ps] [DAP dv (2.104) 


It may be recognized that [DH] is the strain-displacement relation given by equation 
2.74 which corresponds to the linear contribution of the stiffness matrix given in 
equation 2.97, and the contribution of the non-linear strains results in the 2nd 


Piola-Kirchhoff matrix, [s], as 
(s] = [P's] [Aza]” [E] (2 [Ara] + [Aco)) {Ps} (2.105) 
[I's] is given in equation 2.79. It may be shown that the matrix [s] takes the form, 
[si] [0] [0] 
[s] =| [0} [s:] {0} (2.106) 

0} {0} [si] 


The expression for the second term in the right hand side of equation 


2.63 is evaluated in the same manner as the linear and nonlinear parts. On using 


equations 2.100 and 2.101, we obtain, 
[ A8e}7Cis}Pau = f {4} (2 [Bus]? + [Brol”){E] (2[Bra}+(Br0)] {4} °dv (2.107) 
It may be seen that by defining 
{“+99F} = f (2[Bul” + (Bool) [E}(2{Brs] + [Beo)] {a}°dv (2.108) 


the expression {d}7 {torr } represents the work done by the external loads at time 
t+ At. Noting that 
{6't4'e} = {5e} (2.109) 
We approximate for the second term in the right hand side of equation 
2.62 such that 
i, _{6e}7{'s}°do © [saretegt (ats)dy (2.110) 
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which represents the internal virtual work, so that the right hand side of the equi- 
librium equation 2.63 is the difference between the external and internal virtual 
work. On substitution of relations 2.97, 2.102, 2.103, 2.107, and 2.109 into 2.63 
and applying the principle of virtual displacement as shown previously, we arrive 
at the incremental equilibrium equation 2.96, which may be solved by Newton-type 
methods. [Ford and Stieman, 1987] 


E. STRESS-STRAIN RELATIONS 
In this section, the stress-strain relations for a composite material, to be used in 
the three dimensional analysis is developed. Two approaches, one based on classical 
laminate theory and the other based on anisotropic material constitutive relations, 
are presented. 
1. Classical and Higher Order Laminate Theories 
Typical structures composed of composite materials are built using sev- 
eral number of laminae, forming a laminate. Each lamina consists of, typically, 
uniaxial fibers embedded in a matrix, such as epoxy-resin, forming a thin plate. 
Figure 2.4 shows principal material axes, labelled 1 and 2 in directions parallel to 
and normal to the fibers, respectively. It may be noted that in each lamina, there 
exists a state of plane stress, as shown in Figure 2.4. 
Assuming elastic orthotropic material, (i.e., the lamina possesses a plane 


of elastic symmetry parallel to the x-y plane), the generalized Hooke’s law may be 


written as 
O71 Qu Q12 Qs 0 0 0 €1 
O2 Q2 Q23 0 0 0 €2 
oy | _ Qx 0 0 0 €3 
oe a ee : (2.111) 
Os symm 2Qss 0 és 
6 2Q66 € 
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Figure 2.4: Lamina coordinate system (2-D) 
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where the plane-stress elastic constants are given by: 


E, 
Qu = —— 
1 — 2 v21 
be WY2E2 
2 «= — 
1 — 2 Vai 
E, 

Q2 = Q3 = ——— 
1-2 va 


Qaa = Gos, Qss = Gis, Qes = Gia 


Q13 = Q23 =O 


(2.112) 


The subscripts 1, 2, and 3 correspond to normal stresses or normal strains 


while 4, 5, and 6 correspond to shear stresses or tensorial shearing strains in yz, zz, 


and zy planes, respectively. 


The stresses in the material coordinate axes are transformed to reference 


coordinate axes (z, y, z) by the following equation: 


Or m? n? 0 
Oy n? m? 0 
0; _ 0 0 1 
Oy: { | 0 0 0 
Oz 0 0 0 
Ory mn —mn 0 


where the direction cosines m and n are given by m = cos@ and n = sin 8. 


0 


0 
0 
n 
m 
0 


—2mn 
2mn 


2 


0 
0 
0 


(m? — n?) 


6 


(2.113) 


The strains may be transformed in a similar manner. Introducing the 


strain transformation, along with equation 2.110 into equation 2.112 results in 


Or Qu Qi2 0 
Oy Qi Q2 0 
Or —_ 0 0 Q33 
Oys “1 0 O 0 
Cis o 0 0 
Ory k Qis Qo 0 


where, 


© 
i 


0 0 

0 0 

0 0 
2Q44 2Qas 
2Q4s 2055 

0 0 
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2Q 16 
2026 


y 
2066 


Qiuim‘ + 2(Qi2 + 2Q66) m?n? + Qaon* 


(2.114) 








Qi = (Qua + Qr — 4Qee) m?n? + Qua (m4 + n'‘) 

Q22 = Qun* +2(Qi2 + 2Qc6) m’n? + Qrm* 

Qs = Qs 

Qis = Qumn — Qumn® — (Qi + 2Qs6) mn (m? — n?) 
Qrs = Qumn® — Qamn + (Qi2 + 2Qe6) mn (m? — n”) 
Qee = (Qtr + Qu — 2Q12) m?n? + Qeg (m? — n?)’ 

Qea = Qaamn? + Qssn? 

Qas = (Qss — Qaa) mn 


Qss = Qaan? + Qssm? (2.115) 


The stress strain relations presented correspond to k‘* lamina. Now, 
consider a laminate composed of N laminae for which, each lamina has a different 
orientation (@), with respect to the laminate z and y axes. For linear elastic plates, 


the function2! form of the displacement may be assumed to be 


u(z, y, z) uo(z,¥) se zu,(z, y) 
v(z, y; z) = vo(z,y) + zv;(z,y) 


w(z, y) = wo(z,y) (2.116) 


and, the linear strains are given by, 


1 
65 = gui + U;,) (2.117) 
so that, 
€xr = Uo,r + ZUi 2 
Cyy = Voy +201 y 
€z: = 0 
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Gz = 5 (v1 + Woy) 


l 
Crz = g(t + Woz) 


1 z 
€xy = 9 (Uo + Vo,r) + 5 (tw + V1,r) (2.118) 


where ug, Uo and wo are the midplane dispiacements, u; and v; are related to the 


rotations of the normals. It may be noted that the in-plane strains, 


1 
€ro = Uo,r, €yo = Voy, Cry = 5 (Uow + vo,2) (2.119) 
and the curvatures are given by 
1 
Ky = Uns, My = Uy, Kry = 5 (tw + U1,2) (2.120) 


We define stress resultants for plate/shell type structures in terms, of 


stresses and shears (see Figure 2.5) as follows for the k‘* layer: 


z 
No = [2a 
& 
Q: = Oyz2dz 
ie , 
A 
M, = [i o=2de (2.121) 
a3. 


Similar expressions are applicable for Ny, Nzy, Q,, M, and M,,, where A is the 
lamina thickness. 


By summing all laminae over the laminate thickness in the following man- 


tae 


ry 


ner, 


€z9 k . Kr 
je a.) = jas se om ry st (2.122) 


which, in matrix form, may be written as 


{N} = [A] {eo} + [B] {x} (2.123) 
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Figure 2.5: Stress resultants and couples in a lamina 
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where, 


rf 
ue 


(Qi;), (he — hes) 


k=1 
N 
= > (i), (42 - 42.) (2.124) 
k=1 
with 7,7 = 1, 2, and 6 and the moment resultants are given by 
{M} = [B] {eo} + [D] {x} (2.125) 
where, 
i 4 
Di = 3 yi (Q.;), (42 - AB.) 2.126) 


with 7,7 = 1, 2, and 6. 

The displacement field, as stated in equation 2.105 is linear in the thick- 
ness direction, resulting in constant shear stresses. To get better accuracy, a higher 
order displacement field may be used (Reddy, 1984]. 

In order to account for the accurate shear distribution, shape factors are 
used in computing shear energies. These factors are typically obtained by equating 
the shear energies. The procedure is outlined for linear and cubic variation of 


displacement fields. The shear energy due to transverse shear stresses is given by, 


1 h 
U.=5 /, [ A (zz (2€r2) + Sys (2€yz)) dz dA (2.127) 


where A is the area bounded by the lamina surface dy. On using Hooke’s Law, we 
get 
U,= 5 a se Gres (2ere)? + Gys (2tye)*) dz dA (2.128) 


Equating this to the linear displacement field and simplifying, 


U, = (5 V [Ges (ur + Woz)’ + Gyz (v1 + woy)”| dA) h (2.129) 
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Introducing a higher order displacement field yields a more realistic stress distribu- 


tion, but in doing so, a shape factor is introduced to yield consistent shear energy. 
Assuming a displacement field in which the displacements are expanded 
as cubic functions of the thickness coordinate, while the transverse displacement is 


assumed to be constant through the thickness, yields 


Uley2) = Uo(zy) + ZUi(ey) + ZU2(z,y) + 2° Us(z,) 
V(zy.z) = Vo(zy) + 2V1(z,y) + 2? vx(2,y) + 2°u3(z,y) 
Wizyz) = Wo(z,y) (2.130) 


With uo, vo, and wo being the displacements of the midplane, the tensorial 


shearing strains are evaluated as, 


[ux(e.) + 2ugzy)z + 3ug(zy)2” + wo.s| 


2e-2 


2€yz [ore + 2vo2y)z + 3us2y2° + wo] (2.131) 


Using the condition that the transverse shear stresses vanish on the plate 


top and bottom surfaces, we have, 


h 
Or: (z, y; +3] = Oy (=, y, 5) =0 (2.132) 
or, 
h h 
€xz:(Z,y, tr] = &|7z, y, tz] =0 (2.133) 
2 2 
Substituting relations 2.132 into 2.130, we obtain, 
uz = w= 0 
uz = ~ 3R2 (wo,2 + U1) 
4 
v3 = —5y2 (woy + 1) (2.134) 
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The displacement field then becomes 


4 
u = ugptzu,— 2a (wor + U1) 

4 
v = utz- ire (woy + 11) (2.135) 
w = Wo (2.136) 


and the shearing strains are given by 


€z2 = 5 (us + woz) f —4 (3) | 


za ya 4(2), 2.137 
fy. = 9 (M1 + Woe SG (2.137) 
Yielding 
1 2 2 + z\? , 
Uo (Bf, [aor + wan? + Gute +s 44) fi ft (GY ] a 
(2.138) 


The shear energy ratio of the two displacement field is found to be 3, so that the 
correction factor for constant shear stress using cubic displacements is 2. 

It is clear from the discussion that the classical plate theory stiffens the 
plate by not taking into account the higher order terms. If we use the higher 
order theory, we need to introduce a correction factor to the shearing strains of the 
magnitude J. Using equation 2.113, the transverse shearing stresses for the k** 


layer are given by 
Orzz, = 2055, €22 ns 2045, €yz 
Oyz, _ 2045, €xz + 20) 44, Ey: (2.139) 


and the resultants are obtained using equation 2.122 as 


Q: =2 (Ass€zz tf Ags€yz) 
Qy = 2 (Ags€zz + Agatyz) (2.140) 
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Note, that equations 2.139 and 2.140 are applicable for any displacement field. 


Hence, for the higher order theory presented, 


Ai = = fey f i (4) [1-43 ]4 (2.141) 
Sin 


Bs is), [ae — fh h2 — A? 2.142 
= 2 (2), [be ne sna ( »)| ee) 


or, 


with 2,7 = 4, 5. 
In the present three dimensional solid element, for which only three translational 
degrees of freedom per node are defined, resultants are divided by the corresponding 


thicknesses to obtain the stress-strain relations with 


Ay Aiz2 0 O 0 2Ai6 





Aiz2 An 0 0 0 2Adr. . 
0 0 0 0 0 0 
(z] 0 0 0 24m 2A 0 oy) 
0 O O 2Ag, 2As5 0 
Ais Az 0 0 0 2Age 
For the special case of isotropic material, the material stiffness matrix is given by 
A+2G r rN 00 0 
r A+2G A 00 0 
r r A+2G 0 0 0 
[E] 0 0 0 G00 (2.144) 
0 0 0 0G 0 
0 0 0 00G 
where 
a (2.145) . 


~ (1+) (1 — 2v) 
2. Three-dimensional Anisotropic Theory 
As an alternative to using laminate theories to obtain [E] matrix, we may 


use anisotropic definition of the laminates. 


a = (Qi) 6 (2.146) 
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These relations are approximated by obtaining the Q,; in the laminate x-y axes by 
suitable transformations and transverse properties (thickness direction) correspond 


to the matrix characteristics. 


F. CONSISTENT LOADS 

In this section, we consider the element nodal loads vector, due to applied 
loads. Using the virtual work principle, the distributed loads, such as surface loads 
and body forces, are converted into discrete loads applied at the element nodal 
points. Discretizing the distributed loads along these lines are referred to as con- 
sistent or work-equivalent loads. Consider a case where a uniform distributed load 
acts on a prescribed face of the element, as seen in Figure 2.6. The consistent load 


vector may be written as 
{r,} = [wy {f*}ds (2.147) 


For uniform distributed surface loads, we have 


{f?} =p {1} (2.148) 


It is worth noting that the interpolation functions on a given surface, say t = 1 
reduces to that of a plane biquadratic Lagrangian ‘scparametric element and are 
presented in Table 2.1. Invoking symmetry, we observe that the forces at nodes 1, 
3, 5, and 7 are equal, and similarly, forces at nodes 2, 4, 6, and 8 are equal. On 


using the shape functions for the t = 1 surface, we have, 


rp = Lf. (1 -r?) ( Salata 
ail 4 
m4 = [ ie 5 (1 - 1?) (1—s)p dr ds — 579 = ap 
1] 1 1 l 
rT) = de fi 1+r)(1-—s)pdrds —~ 98 — 572 — Gls = oP (2.149) 
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TABLE 2.1: SHAPE FUNCTIONS FOR 9 NODED BIQUADRATIC 
ELEMENT 





N, = 7 (1+ r) (1-5) ~5Ne— 5M ~ 2No 

N3 = z(ltr)(1+s)— 5M — 5M =No 

Ne Ss i(1—r) (+9) ~ 5Na— 5Ne— To 

N, = j(1—F) (1s) — 5a 3Ne — TNo 

N, = 5 (+r) (1 - $) — SNe . 
Ny = 5 (1- r?) (1+8)—5No 

Nee = (1- r)(1- 3) — 5 No 

1 


Ne = =(1- r?) (1-5) - 50 
M = (1-¥) (1-2) 


Note: Node numbering is referred to Figure 2.2 where ¢ = 1, upper plane. 


as NO 
_ 


Figure 2.6 gives consistent element nodal loads for a single element. As a check, the 
total pressure loading on the surface, 2 x 2 x p = 4p, is seen to be equal to the 
sum of all the discz2t'zed nodal point forces. The procedure may be extended for 
more than one element by summing loads at joint nodes, as illustrated in Figure 2.7 


for four elements. 


G. INTEGRATION 
In this section, we summarize the Gauss method for numerical integration, 


including a discussion on some aspects of integration schemes. 
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Figure 2.6: Consistent loads 
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Consistent loads 


Fig: 2:6; 











Figure 2.7: Consistent loads on a mesh 
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Consistent loads on a mesh 


Fag. 2.7% 








1. Gauss Quadrature 
The nature of finite element matrices suggests the usage of numerical 
quadrature. Gauss’s integration scheme is the most commonly used approach and 
is adopted in the present analysis. 
The method enables exact evaluation integrals, consisting of polynomi- 
als of any order, by using appropriate order of integration. In general, the Gauss 


quadrature for a function ¢(r,s,t), has the form 


1opl pl 
I= i, 2 I ¢(r,s,t) dr ds dt DD Limiwjnd(rs s,t) (2.150) 


The integration limit reflects the limits of non-dimensional ‘master’ isoparametric 
elements, while ¢(r,s,t) represents the stiffness contribution. 

Figure 2.8 demonstrates the application of the method for a two dimen- 
sional biquadratic element. Using the weighting factors as given in Table 2.3, the 
element stiffness matrix is evaluated, for example, by using a 3rd order integration 


scheme as follows: 


58 


[K] = ($1 + $3 + $7 + bg) + 99 (G2 + 4 + 6 + bs) + x5 os (2.151) 


wolo 
wolo 


where 


$; = h[B(r, s)]” [E][B(r, s)} | J(r, s) | (2.152) 


as is evaluated at Gauss point 2 as shown in the Figure. 
2. Integration Scheme 
The term “full integration” refers to an integration scheme which evalu- 
ates the integral exactly as shown in the previous example. In the same manner, a 
lower order integration is referred to as ‘reduced integration’. 
In the present analysis, ‘full integration’ is used to evaluate the stiffness 


matrices. When a crude mesh is used, a stiffer structure is obtained. In general, there 
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Figure 2.8: Gauss points 
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TABLE 2.2; SAMPLING POINTS AND WEIGHTS FOR GAUSS 
QUADRATURE OVER THE INTERVAL -1 to 1 


Order n Location of Sampling Point Weight Factor W; 
1 0. 2: 
2 +0.57735 02691 89626 = +i l, 


3 +0.77459 66692 41483 = +/0.6 0.55555 55555 55555 = 5 


0. 0.88888 88888 88888 = § 
4 +0.86113 63115 94053 = + [242]? 0.34785 48451 37454 = }-— + 


+0.33998 10435 84856 = + [357]? 0.65214 51548 62546 = 1+ 2 


r 


where r = V1.2 and (2r — 1) is the polinom order 


are two ways to soften the structure. One way is to refine the mesh and another by 
using ‘reduced integration’. Thus, by using a ‘reduced integration’ scheme, a faster 
convergence and more cost-effective, accurate solution may be obtained. However, 
the method suffers such drawbacks as mesh instabilities or mechanisms, resulting in 


a singular element stiffness matrix. 


H. BUCKLING ANALYSIS 
1. Introduction 
It is well known that thin columns or plates under axial compression tend 
to buckle. Elastic buckling occurs when the compressive stress is well below the 
material stress limit. A flat plate under axial compression shortens in the direction 
of the applied compressive loads. This shortening results in coupling between in- 


plane and out-of-plane displacements. 
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As the applied compressive load increases, there is a configuration at 


which the plate offers no more resistance to deform, resulting in a state of neutral 
stability. The load corresponding to this configuration is referred to as the buckling 
load and constitutes a limit point on the load response curve. 

At this critical value, the deflection becomes very sensitive to any change 
in the configuration. For some structures, beyond the limit point, the load-displace- 
ment path may take any of multiple paths. The point where the plate can take any 
of the different paths is called the Bifurcation point and is illustrated in Figure 2.9. 
In analyzing for nonlinear response, the incremental load method is adopted, which 
may be summarized as follows: (a) the tangent stiffness matrix is formed, and solved 
for displacements for an incremental load. Keeping the stiffness matrix constant, 
corrections to the incremental displacements are obtained in an iterative manner 
until equilibrium is achieved, (b) total displacements for this load are obtained, 
(c) a new tangent stiffness matrix is formed at this new equilibrium position and 
steps (a) and (b) are repeated. This procedure is continued until the desired load 
is reached or the critical buckling load is reached. 

2. Implementation 

In this section, the Finite-Element formulation for buckling will be pre- 
sented (Bathe, (1982), Kolar, et al., (1985)]. The problem of instability can be 
approached either by looking at the equilibrium of the structure in the deflected 
position and transforming all quantities to the initial configuration or by solving the 
system in the current configuration. The former approach, described earlier as the 
total Lagrangian formulation, is adopted here. By performing an incremental load 


analysis, using the nonlinear formulation described earlier, we may write 


([*Ki] + [‘Kwe]) (ay) =O" a P} - (FY (2.153) 
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Figure 2.9: Instability and bifurcation points 








where {P} represents the total load applied and ‘+4*) is a scalar, referred to as 
load parameter. The value \ scales the incremental load and may be treated as a 
constant or variable during iterations. Buckling load is reached when displacements 
become large with no increase in the incremental load, i.e., the global stiffness of the 
structure, as illustrated for a single degree of freedom system in Figure 2.9, becomes 
small and [K] tend to be singular. Thus, using Newton-Raphson and modified 
Newton-Raphson methods, convergence difficulties are encountered as buckling load 
is approached. This is overcome by using arc length methods, described in the next 
section, where the load parameter is continuously updated to reflect the state of the 
structure. 
3. Constant Arc Length Method [Kolar and Kamel (1985)] 

When using the Newton type iteration schemes, the stiffness matrix ‘be- 
comes singular as limit points are approached. In order to obtain post-buckling 
response, a method to overcome this singularity is needed. This is accomplished by 
treating the load parameter as a variable and thus have an adaptive load incremen- 
tation. This approach differs from the conventional Newton type schemes where the 
load level is held constant for all iterations at a given load step. Symbolically, at 


load step m and iteration i, equation 2.153 can be rewritten as follows 
™ [K] {d}*) = (™A + AM + Ad) {P} - {F}O (2.154) 
where ™ [K] is the danwent stiffness matrix at load step m 
Ath = MO 4 AD (2.155) 


The Arc Length Method (ALM) may easily be visualized for a single degree of 


freedom as shown in Figure 2.10. The displacements are updated as 
{z}"F) =™ {2} + {u} + Au (2.156) 
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such that z+") corresponds to the displacement at the (i + 1)" iteration of load 
step m. 
In the constant ALM, the radius of the arc at each load step is constant. 
It is clear from Figure 2.10 that the ALM is used in conjunction with the modified 
Newton-Raphson method, and may also be implemented with NR iteration schemes. 
The method allows one to obtain postbuckling response but bifurcation problems 
require modification that will seek out multiple paths after a limit point. 
4. Convergence Criterion 
For a given load step, the iterations on displacements are carried out until 
a pre-set convergence is achieved. There are three convergence tests most commonly 
used, (a) Displacement Convergence, (b) Residual Force Convergence, and (c) Strain 
Energy Convergence. These criterion may be summarized as follows: 
{Au'}7{Au'} 
{Aut}7{Au'} 
{9'}7{9} ec 
(AM)? {P}T{P} 


{Au'}"{g'} 
AM {Au }{P} 


QDISP 


IA 


ORF. 


IA 


QDISPaR.F. (2.157) 


It may be noted that {Au'} is the incremental displacement at i** iter- 
ation, {g'} is the residual force at i** iteration, AA!{ P} is the incremental load at 
the first iteration, {Au'} is the incremental displacement at the first iteration, and 
a’s are the prescribed convergence parameters, usually in the order of 107? to 1074. 

It is further noted that the initial load {P} used to start the analysis is 
set arbitrarily and only the load parameter is modified automatically to go from 


zero load to the desired load level. 
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Figure 2.10: Constant Arc Length Method 
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Ill. PROGRAM IMPLEMENTATION 


A. INTRODUCTION 
This chapter presents certain aspects of computer implementation of the ele- 
ments matrices developed in the previous chapter. As mentioned in Equation 2.25, 


the general problem to be solved is given by 


[K] {d} 
[K] {D} 


{r} (3.1) 
{R} (3.2) 


where equations 3.1 and 3.2 are the static equilibrium equations for the element 
and structural assemblage respectively. If the stiffness matrix, [A], is independent 
of displacements, the analysis reduces to solving a set of linear algebraic equations. 
In the case of the stiffness matrix being displacement independent, the structural 
behavior is nonlinear, and an incremental load analysis together with a suitable 


iterative methods has to be adopted. 


B. LINEAR ANALYSIS 

In the case of linear analysis, we assume small displacements and small strains, 
and the resulting force-displacement relations are solved only once. Using the dis- 
placement independent neck in equation 2.61 to get the strain-displacement rela- 
tions [Bzo], as given by equation 2.80, equation 2.20 is used to form the element 
stiffness matrix. A series of Fortran subroutines was developed incorporating the 
element stiffness matrix for this element. The material characteristics may be either 
isotropic, laminate theory definitions, or anisotropic description. The subroutines 


are implemented in an existing computer program, FEMCOM, which is capable of 
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element assembly and subsequently calculation of the displacement solution for both 


linear and incremental load methods. 


C. NONLINEAR ANALYSIS 

In this research effort, geometric nonlinearities, namely, large displacements 
but small strains are considered. Consequently, the element consists of a linear 
displacement independent stiffness matrix, [Ko], and two other contributions. The 
first contribution is due to the linear displacement dependent stiffness matrix, [Kz], 
based on [Bz,], as given in equation 2.81. The other contribution comes from 
nonlinear stiffness matrix, [Kz], as given by equation 2.95. 

Note that the stress-strain matrix, [E], may be used both for isotropic as well 


as composite materials using relations 2.131 and 2.132. 


D. SOLUTION PROCEDURE 
1. Composite Material 
In order to generalize the procedure of implementing the solid element 
with composite materials, the plate built of solid elements may be stacked in all three 
directions. Figure 3.1 shows such a stack, where rows of elements are arranged in the 
thickness direction. For each finite element, the stress-strain matrix is computed in a 
subroutine separately, though it would be more efficient to compute it for the whole 
row of elements, taking into account the appropriate layers. In assigning a certain 
number of layers in each row, a constraint to be noted is that the total number of 
layers of all rows match the number of layers of the structure being modeled. 
2. Linear Case 
As mentioned earlier, the matrices corresponding to the linear displace- 


ment independent part was coded into several subroutines and implemented into a 
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Figure 3.1: Thick composite plate-element arrangement 











general purpose finite element program, FEMCOM. The program does automatic 


element assembly and yields solutions to prescribed loads. The flow chart shown in 


Figure 3.2 shows various steps that may be summarized as follows. 


. The material properties, model geometry, applied loads, integration scheme, 
boundary conditions and solution parameters are input. The material proper- 
ties needed for isotropic material are Young’s modulus and Poisson ratio. For 
composites, data needed includes the number of layers, rows of elements, fiber 
orientations, Young’s moduli, shear modulus in three directions, and Poisson 


ratio. 


. Using the shape function derivatives, the coordinates transformation relations, 


Jacobian and the strain displacement relations are established. 


. Using the specified Gauss quadrature, the element stiffness matrix is formed 


in global coordinates. 
. The element global stiffness matrix is assembled. 
. Using Gauss elimination technique, the displacement vector is computed. 


. Stresses may be computed using equations 2.14 and 2.15. 


3. Nonlinear Case 


In order to obtain nonlinear response, either for studying the extension- 


twist-flexure coupling or nonlinear buckling and post-buckling, the analysis pro- 


cedure is termed the incremental load method, and a variation of Newton-type 


iteration is used. The element formulation, assembly and equation solving proceed 


as before, except that additional element stiffness contributions have to be taken 


into account. The assembly and solution to get displacements needs to be done as 
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FIG. 3.2: Flow chart- 
linear analysis 


Figure 3.2: Flow chart — linear analysis 
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frequently as the load steps increments and iterations continue, depending on the 
solution strategy selected. A typical flow chart is given in Figure 3.3. 

For a given load step, the incremental displacements are computed it- 
eratively until the convergence criterion is satisfied. At that point, equilibrium is 
achieved and new incremental load is applied and a new tangent stiffness matrix is 
computed. The iterations continue until the new equilibrium position is obtained. 
In the modified Newton-Raphson method, the tangent stiffness matrix is kept con- 
stant for all iterations for a given load step, while, for the Newton-Raphson method, 
the stiffness matrix for the whole structure is formed at every iteration. By tracing 
the load-displacement path, critical points, characterizing buckling, and stable and 


unstable regions of post-buckling equilibrium states may be identified. 
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Figure 3.3: Flow chart — Nonlinear analysis 
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IV. NUMERICAL EXAMPLES 


A. INTRODUCTION AND NOTATIONS 
In this chapter, selected numerical examples are used to evaluate element 
idiosynchrasies and demonstrate its application in solving critical structural com- 
ponents that use thick composites. Solutions obtained here are compared with 
available elasticity solutions or other numberical solutions. 
1. Material Properties 
In all the examples to be discussed, the material characteristics used are 
as follows. For isotropic materials, 
E = 30 x 10° psi, v = 0.30 
and for composite materials used for laminated plates, layer properties are given by 
E, = 40x 10° psi 
E, = 10° psi 
Giz = Gi3 = 0.6 x 10° psi 
G23, = 0.5 x 10®psi 
y = 0.25 
An eight-layered symmetric laminate configuration using this material is selected. 
All the dimensions presented in this chapter are in inches. In the discussion on 


effects of numerical integration rules, L x M x N notation refers to the number of 


integration points in z, y, and z directions respectively. 


B. COLUMNS AND BARS 


Two simple cases have been selected as part of the element validation process. 
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1. Bars 
A bar clamped at one end and loaded at the other end with uniformly 
distributed traction (Figure 4.1a) was studied and compared to the theory of elas- 
ticity. In the numerical solution, work-equivalent loads were used. The dimensions 
of the bar are 10 in. x 1 in. x 1 in., and it is isotropic. The boundary conditions 
are given by 
u (0. +5) =v (0. i +5) =w (0,345) =0 (4.1) 
Figure 4.2 depicts the effects of reduced integration and mesh refinement on the 
maximum deflection. It is obvious that when the mesh is refined in the thickness 
direction, for instance, one element in each of z and y direction and two in z di- 
rection {1 x 1 x 2] mesh, provides a stiffer solution than for the [1 x 1 x 1] mesh. 
It may be noted that the full (Ff) and reduced (A) integration schemes converge 
to about 95% of the classical solution (See Appendix C). It may be noted that the 
classical elasticity solution does not account for transverse shear stresses. A reduced 
integration in the axial direction (2 x 3 x 3) gives the same results as the full inte- 
gration. However, when reduced integration in the thickness directions is performed 
(3 x 2 x 2), the solution converges slowly. On using (2 x 2 x 3) integration scheme 
in the thickness direction for (12 x 1 x 1) mesh, spurious mode is observed. Table 
4.1 summarizes the effects of various integration schemes and mesh sizes. 
2. Beams 
The next example considered is a clamped, cantilever beam loaded at the 
free end by a shear load. Using the clamped boundary conditions, dimensions and 
material as the previous example, solutions using full and reduced (R) integration 
are compared with the elasticity solution in Figure 4.3. The comparisons also include 


the solution obtained using eight noded first order solid element of ‘GIFTS’ software. 
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Fig. 4.1; Bar Sample problems 


Figure 4.1: Bar Sample Problems 





Fig. 4.2; Clamped Bar Under Uniaxial Load 
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Figure 4.2: Clamped Bar Under Uniaxial Load 











TABLE 4.1: EFFECTS OF REDUCED INTEGRATION AND RE- 
FINED MESH, ON THE MAXIMUM DEFLECTION OF CLAMPED 
ISOTROPIC CANTILEVER BAR UNDER UNIAXIAL LOAD 


oo Pot d.o.f. Integration 
Unies 
Rule Vetasticity X 100 


a: x 3x3) 97.49 













Micali =Ixlxl 
- x 3x3) 98.17 
(2x2x i 129.56 
=2x1xl 99.43 
ee 
115.98 
eo) as ea 90.77 
91.20 
a 91.83 
3=3x1lxl 165 F 100.40 
R 100.81 
RR 112.64 
4=4x1x 1] 219 F 101.04 
R 101.59 
RR 111.39 






6=6xlxl 327 : 101.83 
102.80 
110.74 
12=12x1lx1 102.74 
105.00 
11 Po RR 1,500.00 _| 00 


Oma AE 


Usiasticity 
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The reduced integration shows a better convergence than both the full integration 
and the first order solid element. It may be noted that the eight noded solid element 
converges more rapidly than the full integration scheme of the present element up 
to about 300 degrees of freedom (d.o-f.). 

It can be seen from Table 4.2 that mesh refinement in the thickness di- 
rection results in reduced performance and one element in the thickness direction 
consistently yields good results. 

In Figure 4.4, the effect of transverse shear deformation is studied for a 
12 x 1 x 1 mesh using reduced integration scheme, and compared to the theory of 
elasticity solution (See Appendix C). The results are summarized in Table 4.3 versus 
the aspect (i) ratio. 

It is clear that for thin beams where the elasticity solution is adequate, 
the present element gives stiff solutions, whereas for thick beams (4 < 10), better 
solutions are predicted. The reason for these effects may be attributed to the trans- 
verse shear stresses. In the case of thin bars, or beams, the element aspect ratio 
is very large and the parasitic shear strains appear at Gauss points, resulting in a 
phenomena called ‘shear locking’ [Cook, 1989, and Hughes, 1978]. When the beam 
is thick and the aspect ratio is of the order of 1/10, the transverse shear stresses 
start to become significant, whereas in the elasticity solution, they are taken into 
account only to a limited degree together with the restrictions of Saint-Venant’s 


principle. 


C. CLAMPED PLATES 
An isotropic clamped plate of dimensions 20 in. x 20 in. x 1 in. under 
a central concentrated load is shown in Figure 4.5a. This problem is studied for 


mesh sensitivity and the effects of different integration rules. The present solution 
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Fig. 4.3; End Loaded Beam Bending 
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Figure 4.3: End Loaded Beam Bending 
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TABLE 4.2: EFFECTS OF REDUCED INTEGRATION AND 
MESH CONFIGURATION ON THE MAXIMUM DEFLECTION OF 
CLAMPED ISOTROPIC CANTILEVER BEAM LOADED AT THE 
END 


Mesh * # d.o.f. | WS Wer ee 100 | WS Wer ee 100 
27 a 8 si t 
(27 solid) 


l=I1x1xl 














3=3x1lxl 
4=4x1xl 
9=9xl1lxl 


12=12x1xl 


Welasticity = 101.0 


F: Full integration 


e 3x3x 3 for 27 solid 
e 2x2-x 2 for 8 solid 


R: Reduced integration 
e 2x 3x3 for 27 solid 


{ 8 solid is generated in “GIFTS”. 
* Mesh configuration for 8 solid is twice of 27 solid in each direction, i.e., 2 x 1 x 3 
for 27 solid is 4 x 2 x 6 for 8 solid. 
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TABLE 4.3: CENTER DEFLECTION VS. ASPECT RATIO (+) OF 
AN ISOTROPIC CANTILEVER CLAMPED BEAM LOADED AT ONE 
END 
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See Appendix C. 
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Fig 4.4; End Loaded Beam Deflection 
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Figure 4.4: End Loaded Beam Deflection 








is compared with the elasticity solution (See Appendix C for details on elasticity 


solutions). 
Invoking symmetry conditions of the problem, a quarter of the plate is modeled 


with the following boundary conditions imposed: 


h h h 

U (0. y; +5) v (0. ¥Y, +5) =w (0, ¥; +5) =0 (4.2) 
h 

‘ («, 0, +3] (4.3) 


and the symmetry conditions. 


u (5, y, z) =e (z, 5: z) =0 (4.4) 


The load was taken as one quarter of the total load. Figure 4.6 shows the comparison 
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of a mesh composed of elements arranged in one row of elements (N x N x 1)*vs. 
a mesh of the type (2 x 2 x M), composed of M rows of elements arranged in the 
thickness direction with 2 x 2 elements in each row. Full integration is employed 
in the computations. Table 4.4 summarizes the resultant deflection and mesh sizes. 
It is clear from this and the previous examples that one row of elements in the 
thickness direction is adequate to predict the response of the structures. Figure 
4.7 presents the convergence characteristics of three integration schemes. It may be 
noted that reduced integration (3 x 3 x 2) in the thickness direction yields very close 
results to that of the full integration scheme. Reduced integration produces good 
results by compensating for the estimation of finite element approximation. The 
in-plane reduced integration scheme (3 x 2 x 2) shows divergence in the computed 
response. It may be mentioned that using one element to model quarter plate 
resulted in much higher deflection than expected. This implies that a one element 
model contains spurious modes and a one-element modeling of plate/shell problem 


should be avoided. On examining the convergence plot, with less than 600 d.o.f., the 
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Fig. 4.5; Plate Sample problems 


Figure 4.5: Plate Sample Problems 
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TABLE 4.4: MESH COMPARISON OF AN ALL EDGES CLAMPED 
RECTANGULAR ISOTROPIC PLATE UNDER CENTRAL LOAD (!P 
= 1000 Ib.) 


eee Db 5e 1 13.8529 
4=2x2xl 4.2230 
9=3x3xl 5.0690 
16=4x4x1l 5.6592 
8=2x2x2 4.0133 
12=2x2x3 3.9955 


16=2x2x4 3.9597 


18=3x3x2 5.0647 





evaluated deflection is within 90% of the elasticity solution. Table 4.5 summarizes 


the effects of various integration schemes and mesh sizes. 


D. SIMPLY SUPPORTED PLATES 


Bending of a simply supported rectangular plate under uniformly distributed 
force is presented herein. (See Figure 4.5b) Both isotropic and laminated plates are 


investigated using one quarter of the plate, as discussed previously. The results are 
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Fig. 4.6; Clamped Plate Under Central Load 
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Figure 4.6: Clamped Plate Under Central Load 
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Fig. 4.7; Clamped Plate Integration Rules 
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Figure 4.7: Clamped Plate Integration Rules 
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TABLE 4.5: INTEGRATION RULES COMPARISON OF AN ALL 
EDGES CLAMPED RECTANGULAR ISOTROPIC PLATE UNDER 
CENTRAL LOAD (1P = 1000) 












Mesh # f. | (3x 3x3) ](3x3x 2) | (2x2x 3) 


d.o. 
145 
325 
577 









16=4x4xl 
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compared with Classical Plate Theory (CPT) as given in Appendix C and higher 
order shear-deformation (HSDT) plate theories [Reddy, 1985, Lo et al., 1977]. 


The following boundary conditions are imposed, 


h 
w (0, Y; +5] = w (2, 0, 4) =0 (4.5) 


and the symmetry conditions are as given in equation 4.4. Reduced integration 
(3 x 3 x 2) is adapted throughout all the computations presented in this section. 
The convergence characteristics of the element and comparison to CPT is depicted in 
Figures 4.8 and 4.9 for both isotropic and laminated plates. In the isotropic case, the 
present element shows convergence within 90% of the elasticity solution for less than 
200 d.o.f. In the case of laminated plates, (Figure 4.9), the classical solution [Vinson, 
1987] gives a more flexible solution than the present element. Table 4.6 summarizes 
the deflections of the isotropic and laminated plates and mesh sizes. It may be 
noted that the classical solution uses laminate theory, which neglects the transverse 
shear stresses, and hence the contribution of these stresses is not taken into account. 
This assumes more significance for thick plates (: < 10 to 15). Furthermore, when 
the plate stiffness in the thickness direction is significantly lower than its stiffness 
in the in-plane direction and when the shear modulus in the thickness direction 
is significant, the classical laminate theory does not predict the response of the 
structure accurately. In the present example, ¢ = 20 and & = 40, Gi3 ® Gio. 
It may be concluded that using laminate theory for bending of thick plates yields 
a nonconservative estimate of deflections and special attention should be given to 
the stiffness ratio [| and shear modulus ratio [2] in determining the response 
of such plates. In Figures 4.10 and 4.11, the maximum deflection is presented for 
different aspect ratios, (2), of the plate. Both isotropic and laminated plates are 


analyzed and compared to CPT. In addition, the solution of the laminated plate is 
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also compared -9 Higher Order Shear Deformation Theory (Reddy, 1985], as shown 
in Table 4.7. 

As in the beam bending cz .e, shear locking is observed for thin isotropic plates. 
For thick plates, (say, t= 4), the computed deflections become significantly larger 
than predicted by CPT, as expected. 

Examining Figure 4.11, it may be deduced that even the HSDT [Reddy, 1985] 
underpredicts the deflections. For thin laminated plates, the shear locking effect is 
not as significant as observed for isotropic plates. This may be attributed to the 
fact the laminated plate has more flexible transverse material stiffness in bending 
than the coefficients than the isotropic plate, so that shear locking is expected to 


develop only for thin isotropic plates. 
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TABLE 4.6: ALL EDGES SIMPLY SUPPORTED RECTANGULAR 
PLATE, UNDER UNIFORMLY DISTRIBUTED LOAD 


Se cou 
Mesh # d.o.f. 


4=2x2x1] 








9=3x3x1l 
16=4x4x1l 


CPT 


* Neglects G,3 and G23. See Appendix C. 
The uniformly distributed load is taken as the total consistent load over the area of 


the quarter plate. 


The maximum deflection is taken at upper surface. 


3 
» = iH) BF 





TABLE 4.7; CENTER DEFLECTION VS. ASPECT RATIO OF SIM- 
PLY SUPPORTED RECTANGULAR PLATE UNDER UNIFORMLY 
DISTRIBUTED LOAD 





*Reference: Reddy, 1985. 
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Fig. 4.8; Simply Supported Isotropic Plate 
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Figure 4.8: Simply Supported Isotropic Plate 


Fig. 4.9; Simply Supported Laminated Plate 
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Figure 4.9: Simply Supported Laminated Plate 
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Fig. 4.10; Isotopic Plate Deflections vs. Aspect Ratio 


2 
jo) 
* 
% 
az 
= 
% 
i) 
~ 
* 
2 
S 
2 


| —e— Mresent Element 
-e- CVT 


WwW 


10! 


Side/Mickness [afr] 


Figure 4.10: Isotropic Plate Deflections vs. Aspect Ratio 


Fig. 4.01; Laminated Plate Deflections vs. Aspect Ratio 
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Figure 4.11: Laminated Plate Deflection vs. Aspect Ratio 


75 














V. CONCLUSIONS AND SCOPE FOR 
FUTURE RESEARCH 


A. CONCLUSIONS 

This study suggests a three-dimensional higher-order finite element to be in- 
corporated in the analysis of thick plates composed of both isotropic and laminated 
composites. By using a tri-quadratic Lagrangian twenty seven noded solid element, 
no assumptions on transverse shear strains are introduced in the formulation. The 
formulation, based on the principle of virtual work, is presented for both linear 
and nonlinear analysis. The material constitutive relations for linear isotropic and 
composite materials are presented. For composites, both laminate theory and three 
dimensional anisotropic adaptations are described. 

Several numerical examples using linear analysis are given for bars/beams and 
plates using both isotropic and composite materials. Three dimensional anisotropic 
relations are adapted for composites. The results show that the present element is 
effective for analysis of thick beams and plates, but exhibits shear locking for thin 
beam and plates. 

Spurious modes are revealed for single element usage in plate modeling, as is 
the case for some other finite elements. 

Reduced Integration in the thickness direction for beams and plates gives sat- 
isfactory results. An interesting outcome is that one element is sufficient to capture 
transverse deformation for thick laminated structures and mesh refinement in the 


other two directions yields convergent solutions. 
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B. SCOPE FOR FUTURE RESEARCH 

More numerical experiments need to be performed to compare the present 
soiution to closed-form solutions [Pagano, 1969] to evaluate the efficacy of this ele- 
ment. 

Implementation of buckling analysis using the nonlinear element matrices pre- 
sented herein is another task that may prove useful in predicting buckling response 
of thick composite cylinders subject to external pressure. By incorporating the cen- 
trifugal force in the external virtual work done by body forces, this element may be 


used in modeling rotor blades. 
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APPENDIX A 


Shape Functions and Derivatives for Solid Element 


Shape Functions for Solid Element 


Mid-edge nodes: Mid-plane nodes: 
No= $r(1+r)(1—s?)t(1+t)  No= $(1—r?)(1—8?)t(1+t) 
Ne= Fr?) s(14e)t(+t) Nar=—b-r4) (1-84) t-9) 
Ne = —gr(1—r)(1—s?)t (1 +8) Niu = qr (+r) (1— 5%) (12°) 
Ng = —j (l—r?) s(1— 8) t (1+?) Nis = —gr(1—r) (1-8?) (1—t?) 
Nio = —4r (1 +r) 8 (1—s) (1-2?) Ni3 = q(1—r°) 8(1 +6) (1—2’) 
Na= fr(ltr)a(t+s)(1-2)  Nivs—$ (1-1?) (1-9) (1-2) 
Nu = —4r (1—r) 8 (1 +8) (1-22) 
Nige= 4r(l1—r)s(1—s) (1-??) 
Noo = ty (1+r)(1—s?)t (1-2) Center node: 
Naz = —5 (1—r?) 8 (1+8)t(1-t) 
No = qf (1—r) (1— 8?) t(1—#) Nig = (1 — r?) (1 — 8?) (1 — ¢?) 
Nog = $r(l—r?)s(1—s)t(1—t) 


M 


N3 


(1+) (15) (144) 5 (Na+ Nat Mio) — 3 (Nai + Mir + No) — Mie 
(1+ r) (148) (148) — 5 (Na+ Nat Min) ~ 5 (Nii + Mis + No) — 2 Mie 
(1 =r) (148) (148) — 5 (Nat No+ Mia) — 5 (Nis + Nas + No) ~ 2 Mie 
(1 =r) (1= 5) (144) — 5 (No-+ Na+ Mis) — 5 (Mis + Mir + No) — 5 Mis 
(1+r) (1 ~ 6) (Lt) ~ 5 (Nao + Nag + Mio) — 3 (Nia + Naz + Nar) - 5 Mis 
(1+ 7) (146) (1-1) ~ 5 (Nao + Naz + Miz) — (Ni + Mas + Nor) — 2 Mie 


8 
. 1 1 1 
(l-r)(1+4)(1-t)- 3 (N22 + Ng + Nig) - q (Mis + N15 + N27) — g Mie 


GO| = GO] = CO] m= GO] m= GO| = GO| m= CO] = CO} me 


1 1 1 
(1-r) (1 — 8) (1—t) — 5 (Naa + Nae + Mie) — q (Mis + Niz + Naz) - g Mis 
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Shape Function Derivatives for Solid Element - r direction 


Mid-edge nodes: 


No, = 4¢(1+2r) (1-5?) (1+?) 
Nas = —f ret (1 +s) (1+?) 

Ne, = _{; (1 — 2r) (1 -- s?) (1+ 2) 
Ne- = t rst (1—s) (1+?) 

Nior = 7 (1 + 2r) (1 — s) (1 -?t?) 


Ni2e = 48 (1+ 2r) (1 +8) (1-?#?) 
Nia = —ba (1 2r) (1 +8) (1-22) 
Niger = 48 (1 —2r) (1-38) (1 —t?) 


Noor = 73 (1 + 2r) (1 — 8?) (1-2) 


Noor = — Great (14+ s)(1-t) 
Nor = ty (1 — 2r) (1— 8?) (1-2) 
Noor = —5rst (1 —s) (1 -t) 


Corner nodes: 


My 
N35 
Nsr 
N75 
Mr 
Nair 
N23,2 


Nose 


Mid-plane nodes: 


Ng,» = —rt (1 — 8”) (1 +t) 


Naz 
Nir 
Nis 
Ni3r 
Nizys 


rt (1 — s”) (1-2) 
1 (14 2r) (1-82) (1- #2) 
f (1 — 2r) (1 — #2) (1 — #2) 


—rs (1+) (1-1?) 


rs (1—s) (1 —t?) 


Center node: 


Niger = —2r(1 — 8?) (1 — t?) 
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1 1 1 
3 (1+s8)(1-t)- 3 (Noo, + Noor + Nizar) — ri (Nie + Nise + Naz) — 


: (1-8) (1+#)- + (Nay + Nae + Nior) - 7 (Mie + Mize + Nor) — 5 Mier 
1 (14s) (140) — 3 (Wap + Nay + Mize) — E (Nite + Mise + Noe) — E Mise 
-; (1 +s) (1+) - 5 (Nay + Nee + Niar) — 7 (Mise + Nis + Nor) — 5 Mise 
-; (l—s)(14t)- 5 (Ney + Nar + Nier) — 7 (Mis, + Nizr + Nor) - 5 Mise 


1 1 . 1 1 
8 (l—s)(1-t)- 2 (Noor + No6r + Mor) — 7 (Nite + Mize + Naz) -— g Mise 


Z 


Mayr 
3 8 


1 1 1 1 
“3 (1+s)(1-t)- 3 (N22,7 + Noae + Niar) — q (Nis + Miss + Noz,r) — =Miar 


8 


1 1 1 1 
“3 (l-s)(1-t)- 3 (Noar + Noor + Nios) — q (Nise + Nite + Nore) - g Mier 





Shape Function Derivatives for Solid Element - s direction 


Mid-edge nodes: Mid-plane nodes: 
No, = —4drat (l+r) (1+) No,, = ~st (1-1?) (1+?) 
Nas = tear) (1+2s)(1+t) Noaz,= st(1—r?)(1-t) 
Noe = 5rst(l1—r) (1+?) Nis = —ra(1+r) (1-t?) 


Ne, =—Lt(1—r2)(1—26) (1+) Nise = rs(1—r) (1-22) 

Nios» =— “Ly (1 4 r) (1 — 28) (1 - #2) Nis, =  4(1—r?) (1 + 2s) (1 -#?) 
Ni20 = te (14 r) (1428) (1-2?) Nize = —£ (L—r?) (1-28) (1- £2) 
—te(1—r) (1428) (1-22) 

Nie = fr(1—r) (1-28) (1-2) 


Noo,= 45rst(l1+r) (1-1?) Center node: 

Noo = —£t (1-2) (1 +28) (1-1) 

Noa, = —4rst (1—r) (1-t) Nia. = —28(1 —r?) (1 —t?) 
Noe, = hear) (1-28) (1-8) 


Corner nodes: 


1 1 1 1 

M,. = =5 (l+r)(1+t)- 3 (No. + Nas + Nios) - q (Mus + Niz,3+ Nos) —- g Mis. 
l 1 1 1 

N3, = 3 (l+r)(1+t)- 3 (No. + Nas + Miz.) — (Mine + Ni3,. + Nos) — g Nie. 


1 1 1 1 
Ngo = g i-r)+8- 5 (Nas + Now + Nias) — 7 (Nise + Nise + Noo) gia. 


Nrs = AE (1) 14) = 5 (Nou + Now + Nig) — 5 (Nise + Nite + Now 1) gM 

Nios = Fosfor Mant ta) List et) Bha 
Nas = ; (1+r)(1-t)- 5 (N20. + No2,5 + M22) - : (Nias + Nis. + N27,2) — 5 Mis. 
N23. = b=) (1=t)— baa + Nana + Mina) — 5 (Nine + Nase + Nara) — 5Nieo 
Nos. = 2 =e) 1-8) — 5 Wane + Nae + Nig) — 5 (Nise + Nita + Nara) ~ 5 Mig 
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Shape Function Derivatives for Solid Element - t direction 


Mid-edge nodes: 
Nar= dt (141) (1—s?) (1420) 
Ner= £a(1—r2) (142) (1 +21) 
Nee = —tr (1—r) (1 — 82) (1 + 2) 
Nee =—£ 4 (1—r2) (1—a) (1 +21) 
Niot = drat (1+r) (1-8) 
Nate —frot (141) (145) 
Nias = grat (l—r) (1+8) 
Nis aired (ese) (Ss) 
Naos = —fr (1-41) (1?) (1-20) 
Naor = — s(1—r?) (143) (1 — 22) 
Nae= ir(1—r) (1-2) (1-21) 
Noe, = ts (1 —r2)(1—s) (1-2) 


Corner nodes: 


Nit 
N3,2 
Ns, 
N74 
M19, 
Nai, 
N23, 


N25, 


Mid-plane nodes: 


No, = 4$(1—r?) (1 — 8?) (1+ 2t) 
Naz. = —4 (1-1?) (1 — 8?) (1 — 2t) 
Nie = —rt (1 +1) (1 — 8?) 


t= —st (1-1?) (1+8) 
st (1—r?) (1-38) 


Center node: 


Nise = —2t(1 - r2) (1 — 6?) 


; (l+r)(1-—8)- + (Na, + Na+ N10) - + (Mins + Ni742+No1)—- 5 Mis 

5 (Ltr) (146) — 5 (Naa + Nae + Nias) - 7 (Mire + Maye + Nox) — 2 Mie 

; (1—r) (1+) - 5 (Nau + Nes + Nias) - 7 (Mise + Nis: + Not) - 5 Nas. 

; (1-r)(1-s)- ; (Nes + Nas + Nie +) — 7 (Mis. + Niza + Not) - aNiss 
-; (l+r)(1-s8)- 5 (Naos + N62 + Ni0,1) — 7 (Mins + Niz,2 + Noz2) - 5 Miss 
=} (144) (149) — 5 (Nao + Naas + Nias) ~ 5 (Nine + Nise + Nar) ~ GMs 
-; (l—r)(1+s8)- 5 (Nan + Noas + Nias) — 5 (Mis + Nis.¢+ No72) — 5 Miss 


1 1 1 1 
“3 (l1-r)(l-s)- 3 (Nae + No6.2 + Nie.) — 1 (Nis,t + N17, + No7,2) — g Mise 
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APPENDIX B 


Jacobian matrix elements: 


Jacobian Matrix 


27 
Jnsat,e, = DL Nip 


t=1 


27 
Ji2 =Yyr = > Nix yi 


t=1 


27 
Jig =2, = > Ni- 2A 


t=1 


2? 
Jy= 25 = > Nis z; 


t=1 


27 
Jn= ye = DNs ¥i 


s=1 


27 
Jn =Zt. = >) Nis ii 


=1 


27 
Jn =2e = DONip 2 


=1 


27 
Jn =yr = DoNit yi 


i=] 


27 
Jg=zs = D Nip Xi 


$21 


Elements of the inverse Jacobian matrix: 


Th 
Ty 
T13 


Pa 


(J22 J33 — Jo3 J32) 
(Siz J32 — Ji2 Js) 


(S12 Jaa — Jia J22) 


Be Pe IR AI 


(J23 J3y = Jo J33) 
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Tas = 


Jacobian matrix determinant: 


J = det [J] 








(Jur Js3 — Jiz Jai) 
(Ja Jiz — Ji J23) 
(Ja1 J32 — Jo2 J31) 
(Ji2 Jai — Jur Js2) 


(Jar J22 a Ja Ji2) 


y 
J 
Z 
J 
t 
J 
Ls 
J 
i 
J 


= Jy (J22 Jaz — Jo3 Js2) 
Jy2 (Jar Ja3 — J23 Ja) 
Ji3 (Ja J32 — Jaz Jn) 





APPENDIX C 


Theories 


A. THEORY OF ELASTICITY SOLUTIONS 
1. Cantilevered bar under traction 


e P = Total load 

e L = Bar length 

e A = Cross section area 

e E = Young modulus : 


2. Cantilevered Beam under end load 


Wmozr = 367+ 976 
PL? 3 h\? 
~ 3EI iegaen($) | 





Reference: Timoshenko, 1951. 


B. CLASSICAL PLATE THEORY (CPT) 
1. All edges clamped rectangular isotropic plate under central load 


Pa? 
Wrnaz = we 
aD 
2 
p= —ER_ 
12(1 ~ v?) 
a = Q.00560 for vy = 0.3 


Reference: Timoshenko, 1959. 
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2. All edges simply-supported, rectangular plate under uniformly 
distributed load 





Wher = er. 


RQ 
It 
ran) 
° 
S 
ra 
ro) 
a 
o 
“t 
y 
Hl 
= 
rae) 


3. Composite 


D = Dym* +2(D2+2De) (rua)? + D2! 





TABLE C-1 
SAMPLE COMPOSITE MATERIAL DATA 


Table C-1; Sample Composite Material Data 


INPUT DATA ; 


LAMINA; THKNES ; THETA ; El ; 





0.40000E+08 ; 
0.40000E+08 ; 
0.40000E+08 ; 
0.40000E+08 ; 
0.40000E+08 ; 
0.40000E+08 ; 
0.40000E+08 ; 
0.40000E+08 ; 


0.10000E£+07 
0.10000E+07 
0.10000E+07 
0.10000E+07 
0.10000E+07 
0.10000E+07 
0.10000E+07 
0.10000E+07 


0.60000E+06 
0.60000E+06 
0.600008+06 
0.60000E+06 
0.60000E+06 4 
0.60000E+06 
0.60000E+06 
0.60000E+06 


OUTPUT DATA ; 


A(i,j)-MATRIX 


0.12773E+08 0.55940E+07-0.81226E+06 
0.55940E+07 0.17603E+08-0.30995E+07 
~0.812266+06-0.30995E+07 0.59436E+07 


B(i,j)-MATRIX 
0.00000E+00-0.27344E-01 0.00000E+00 
-0.27344E-01-0.16406E+00 0.00000E+00 
0.00000E+00 0.00000E+00-0.62500E-01 

D(i,j)-MATRIX 
0.20743E+07 0.28699E+06 0.72461E+05 


0.28699E+06 0.81544E+06 0.17998E+06 
0.72461E+05 0.17998E+06 0.31612E+06 


a a a a a a a a a a a i we a a ww wr we ewe 


Note; A,B and D matrices are evaluated ,neglecting Transverse Shear 


eexe contribution i.e. G13=G23=0 


Using Navier’s solution with n=m=200, 


wmax = 0.052328 
3 
Wmax*E *h 


2 #100 = 0.3634 
we) 06h 


Where q=90 ; a=20 
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i.e. 100 terms for each 
dicection, as given in the above, we have, 








10. 


12. 


13. 
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